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When a shallow layer of inviscid fluid flows over a substrate, the fluid particle tra- 
jectories are, to leading order in the layer thickness, geodesies on the two-dimensional 
curved space of the substrate. Since the two-dimensional geodesic equation is a two 
degree-of-freedom autonomous Hamiltonian system, it can exhibit chaos, depending 
on the shape of the substrate. We find chaotic behaviour for a range of substrates. 
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I. INTRODUCTION 

Many well-known physical systems take the form of geodesic flow on a manifold. For 
instance, Euler's equation can be though of as a geodesic flow in the space of volume- 
preserving diffeomorphisms, [1-6] and free rigid body motion as geodesic flow in SO (3). In 
both cases, the metric on the space corresponds to the kinetic energy norm. It is also known 
that the geodesic deviation equation [7] describes the stability of such flows. For instance, 
a space of negative curvature will lead to divergence of trajectories, and hence to chaos if 
the space is compact. But compact spaces of strictly negative curvature are hard to come 
by in the real world, to say the least. If we expect the negative curvature to lead to chaotic 
geodesies, we are better off looking for spaces with non-sign-definite curvature, but such that 
the averaging of the curvature over trajectories leads to chaos {i.e., the negative curvature 
'wins'). 

In this contribution we will discuss a system which is physically-motivated and leads to 

chaotic geodesies. This system is the flow of a shallow layer of ideal, irrotational fluid on a 
curved substrate. Following Rienstra [8], we will show that, to leading order, the governing 
equation can be solved in terms of characteristics. Moreover, the characteristics are geodesies 
on the curved substrate, possibly modified by gravity if it is present. 

Of course, the chaotic trajectories have a nasty tendency to cross and form caustics every- 
where. Hydrodynamically, caustics are usually manifested as hydraulic jumps or so-called 
'mass tubes,' [9] visible as a thicker edge region of the fiuid in Fig. 1 (top). Edwards et al. [9] 
have recently described these mass tubes using the theory of 'delta-shocks.' [9-13] At this 
point, we are unable to apply this theory to our problem, which means that solutions be- 



*To appear in Proceedings of the conference on Chaos, Complexity, and Transport: Theory and Applications 
(Le Pharo, Marseille, June 2007). 
^Electronic address: jeanluc@math.wisc.edu 



FIG. 1: Experiment in the kitchen sink, using a cut-open plastic bottle (top). The jet from the 
faucet impacts the inclined bottle. The pattern is qualitatively well reproduced by following fluid 
trajectories emanating from a point source on a cylinder (middle). However, if we pursue the 
trajectories further (bottom), the ideal theory presented here predicts that the flow should crawl 
back up the side to the same initial height. The discrepancy is clearly due to dissipation. 

come dubious after characteristics begin to cross. Unfortunately, since our characteristics 
are chaotic, they tend to cross a lot. Nevertheless, we believe that studying the basic prop- 
erties of this geodesic flow is worthwhile as a first stab at describing the transport properties 
of flows on curved substrates. In addition, the geodesic flow we present is an interesting 
mathematical system, with rich dynamics that deserve to be studied on their own. 

Having gotten these disclaimers out of the way, let us proceed with the analysis. We hall 
do this in stages. In Section II we introduce a curved non-orthogonal coordinate system to 
describe the substrate, singling out a direction normal to the substrate. In Section III we 
use this direction to expand our fluid equations and derive a shallow-layer form. We show 
that the resulting equation can be solved in terms of characteristics, which are geodesies on 
the two-dimensional curved substrate, modified by gravity. In Section IV we look at specific 
numerical solutions for particle trajectories, and in Section V we speculate on their chaotic 
nature. We offer some closing comments in Section VI. 
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II. COORDINATE SYSTEM 
A. Separating the Shallow Direction 

In our problem, fluid motion occurs over a curved substrate of arbitrary shape. The 

direction normal to the substrate is special in that it defines the direction in which the fluid 
layer is assumed 'shallow.' Hence, it is convenient to locate a point r in the fluid as 

r{x\x^, y) = X{x\ x^) + y e^{x\ x^) (1) 

where X{x^, x"^) is the location of the substrate, 63 is a unit vector normal to the substrate, 
and y is the perpendicular distance from r to the substrate. The coordinates and x'^ are 
substrate coordinates used to localise points on the substrate. For example, in Section II B 
we will use the Monge parametrisation, X = (x^ x"^ , where / gives the height 

of the substrate. 

The tangent vectors to the substrate are 

ea := daX (2) 

where da '■— d/dx". The coordinate vectors associated with the coordinate system are found 
from (1), 

dv 

ea — daV = ea + 0{y) , 63 := — . (3) 

where Greek indices only take the value 1 or 2. Note that the are not necessarily 
orthogonal or normalised. 

We adopt the convention that quantities with a tilde are evaluated in the 'bulk' (away 
from the substrate), and thus depend on y, whilst those without the tilde are 'substrate' 
quantities and do not depend on y. Thus, ea{x^,x^,0) — ea{x^,x^). The three-dimensional 
metric tensor Qab has components 

0a/3 — ea • 6/3 = Ga/S , 0a3 — • 63 = , ^33 := 63 • 63 = 1 , 

where 

:= e„ ■ eg = G^/s + 0{y) , 

The three-dimensional metric tensor is thus block-diagonal, and the y coordinate is un- 
stretched compared to the Cartesian coordinate system. It measures the true perpendicular 
distance from the substrate to a point in the fluid. 

Given the substrate vectors e^, it is easy to solve for the covectors e", which are such 
that e"^ • = 5/3°. Then the bulk covectors are 

e" = e" + 0(|/), (5) 

to leading order in y. From (5), we find the inverse metrics. 



Higher-order terms in y will not be needed. 
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B. Substrate Coordinates 



For most applications in the literature of thin films and shallow layers, orthonormal co- 
ordinates have been the coordinates of choice. This is because the main substrate shapes 
that have been treated are planes, cylinders, and spheres, where orthonormal coordinates 
are readily available. For a general substrate shape, orthonormal coordinates are difficult 
to construct and require numerical integration. Singularities (umbilics) also cause prob- 
lems. [14] For our application — flow down a curved substrate — the Monge representation of 
a surface [15] is the most convenient. 

The Monge representation is a glorified name for a parametrisation of the substrate by 

X{x\x^) = {x^ x^ f{x\x^)f (7) 

in three-dimensional Cartesian space. Following standard notation, [15] we define 

P = dif, q = d2f, (8a) 
r^didj, s^did2f, t^d2d2f. (8b) 

The unnormalised, nonorthogonal tangents Ci and 62 are 

ei = 9iX = (l Q pf , e2 = 92X=(0 1 qf , 

and their normalised cross product gives the normal to the substrate, 

1 . .xT 



63 = - (-P -q 1) 



The corresponding covectors are 



and 63 is its own covector. The metric tensor of the substrate and its inverse are 

= e„ • e, = +f ^ ^^^2) , (9a) 

{G-^} = e'^ . = i- f ^ + ^\ , (9b) 

\ -pq 1 + p y ^ ^ 

w = (det G„^)'/' = y/l+p^ + q^. 
Finally, we will need the Christoffel symbols F^^, defined by 

= l'^"'' (daG^P + dfsG^c - d.Gap) (10) 
and in Monge coordinates given by 



with determinant 
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The Christoffel symbols arise when taking covariant derivatives. [7, 16, 17] Note that in (10) 
we used the usual convention that repeated indices are summed. 
We write the normalised gravity vector as 

gr=(sin^cos0 sin ^ sin —cos 6')"^, (12) 

so that the inclination angle 9 is zero for gravity pointing downwards, and for G 
(— 7r/2,7r/2) positive 9 induces flow in the positive direction. Then we have the com- 
ponents 

Qs^ = 9 • = — {p cos 9 + pq sin sin — (1 + g^) sin cos 0) /w'^ , ^^^^ 
= g ■ = — {q cos 9 + pq s'm9cos(f) — (1 +p^) sin sin 0) /w'^ , 

The specific parametrisation of the substrate introduced in this section will not be needed 
in the derivation of the equations of motion (Section III), only in their solution. Hence, a 
different parametrisation could be used if called for by the geometry of the substrate. For 
instance, flow down a curved filament is better parametrised by cylindrical coordinates, or 
if the substrate has overhangs (making / multivalued) coordinates based on arc length are 
preferable. 



III. EQUATIONS OF MOTION 

Now that we've set up an appropriate coordinate system on our curved substrate, we need 
some dynamical equations of motion for the fluid. We assume an inviscid, irrotational fluid 
with a free surface at y = ri{x^,x'^), with slip boundary conditions at the substrate y = 0. 
The pressure on the free surface is assumed constant (zero). We also assume the flow 
is steady and irrotational, so that the the velocity can be written in terms of a scalar 
potential, u — V(/?. The equations satisfied by the fiuid are then 

V^V^ = 0, mass conservation; (14a) 

l\V(pf + -^-g -r = H, Bernoulh's law; (14b) 

where H is a constant, with boundary conditions 

dycp — at y = 0, no-throughflow at substrate; (15a) 
V</? • V77 = dy(p ai y — r), kinematic condition at free surface; (15b) 
p — at y — T], constant pressure at free surface. (15c) 

In terms of our curvilinear coordinates, equation (14b) becomes 

G"^ d^^ dpif + {dy^f + ^-2g.r^2H. (16) 

A. Small-parameter Expansion 

Now we assume that the fluid layer is shallow, so that y is proportional to e. After 
replacing y by ey, Eq. (14b) becomes 

(G"'' + 0{e)) dc^ipdpip + e-\dy^f + - 2g ■ {X + ey e^) ^ 2H. (17) 
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FIG. 2: A pencil of 30 trajectories starting from the origin at different angles, each with initial 
kinetic energy 1/2. The substrate has shape f{x^,x^) = /o cos cos x^, and the plots are for 
different values of /q. Gravity is turned off: g = 0. The grey background shows the periodicity of 
the substrate. 



We also expand (p in powers of e, 

(p{x^, y) = (/?(o) + e + ip(2) + ■ ■ ■ ■ (18) 

The leading-order term in (17) occurs at order £~^, and gives 9j/</?(o) = 0. Hence, we 
have </7(o) = $(x^,x^) independent of y. The next nontrivial terms are at order £°, 

G'^^da^ dp^ + {dy^{i)f - 2g ■ X = 2H. (19) 

We evaluate the whole of (19) at y = 0, and use the boundary conditions (15a) and (15c), 

G'^^da^ dp^-2g-X = 2H. (20) 

This is the equation that we need to solve to find the leading-order velocity poten- 
tial ^{x^,x^). We discuss the method of solution in the next section. Note that we will 
not need to solve the mass conservation equation (14a): at leading order, it only serves to 
find the fluid height once the velocity field is obtained. 



B. Solution in Terms of Characteristics 



As pointed out by Rienstra, [8] the trick to solving Eq. (20) is to use the method of 
characteristics. To do this, we differentiate (20) with respect to x^, which gets rid of the 
constant H, 

2G"^9a$ d^dp<^ + d^G^f^da^ d^^ = 2g ■ d^X. (21) 
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The horizontal components of velocity are = G^'^dp^, where the overdot denotes a time 
derivative; hence, 

d^dp<l> = dpiG^si:^) = G^sdpx^ + dpG^sx^ . (22) 

From the chain rule, we have = dpx^x^; using this and (22) in (21), we find, after dividing 
by two, 

G^sx^ + dpG^sif^x^ + Id^G'^'^GasGpp x^x^ ^ g e^. (23) 
Now we multiply by G'^'^, and obtain after an integration by parts and a bit of manipulation 

x" + r''^px"x'^ = g-e'' (24) 

where the F^^ are defined by (11). Equation (24) describes geodesies in the curved coor- 
dinates of the substrate, under the influence of gravity. In the absence of gravity, the fluid 
trajectories are essentially going in straight lines in the curved substrate coordinates. (In 
general relativity, unlike here, the gravity determines the curvature of space.) 

If we define the covariant derivative of a vector along the trajectory [7, 16-18], 

D 



V -.^V" ^Vlf^x'^V^ , (25) 



where r is the time (to avoid confusion with t in Eq. (8)), then the geodesic equation (24) 
takes the more intuitive form 

^x'^ = ge^ (26) 

which looks a lot like Newton's second law, but here it incorporates the constraint that fiuid 
particles remain on the substrate. 

Equation (24) is a two degree-of-freedom autonomous Hamiltonian system, with the 
energy H defined by equation (20) as an invariant. Hence, any other invariant will make the 
system integrable, and rule out chaos. In particular, a surface with a translational symmetry 
cannot exhibit chaos. 



IV. FLUID PARTICLE TRAJECTORIES 

To get a feel for the possible range of behaviour of fluid trajectories, we now solve the 
geodesic equation (24) for a range of substrates. We shall always use a set of fluid trajectories 
starting at the same spatial point, with the same initial kinetic energy but different direction. 
This models a point source, or a thin jet impacting the substrate. 

First, following Rienstra [8] we solve the equations on a cyhndrical substrate. Figure 1 
(middle) shows some trajectories, all emanating from the same point. Qualitatively, the 
pattern captures well the observed behaviour of a jet (from a faucet) impacting the inside of 
a cut-out plastic bottle (Fig. 1, top). However, if we pursue the trajectories further (Fig. 1, 
bottom) , we see that they crawl back up the side of the cyhnder, with no loss of energy, in 
contrast to the experimental picture. This comes from neglecting the hydraulic jumps that 
occur, as well as viscosity. [9] Observe that the trajectories follow a very ordered pattern, 
and are definitely not chaotic. This is as expected, since there is a symmetry direction, and 
so the motion is integrable (Section HI). 

Next we move on to more complex substrates. Since there is basically an infinity of 
choices here, we limit ourselves to periodic substrates with shape 



/(x^,x^) = /o sin x"*^ sin 



(27) 
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FIG. 3: Same parameters as in Fig. 2, but with gravity turned on: g = 1- The trajectories exhibit 
chaos-hke behaviour for much lower substrate height, since they begin at the top of a bump and 
thus have potential energy to draw upon. 



for various values of /q. The other variables in the system are the strength of gravity (which 
can be chosen as unity if it is not zero, by rescaling the substrate height) and its orientation 
(as given by the angles 9 and (f) in Eq. (12)). 

Figure 2 shows a pencil of 30 trajectories starting from the origin at different angles, 
each with initial kinetic energy 1/2, for different values of /o, in the absence of gravity. 
The first two cases display regular behaviour, but for substrate heights /o = 0.7 there is 
chaotic-like behaviour. These are, however, fairly extreme values of /o, corresponding to 
heavily-deformed substrates. Our expansion should be able to accommodate this, since the 
variations in the substrate height are not assumed small (only those in the fluid thickness 
are). For extreme heights (/o = 1.2, last case in Fig. 2), some trajectories actually backfire 
and come around the initial point. 

Figure 3 shows results for the same parameters as Fig. 2, but with gravity (? = 1- The 
inclination is nil {6 = 0). It is clear that chaotic-like behaviour sets in for much smaller 
values of /o, even showing backscatter for /o = 0.5 in the last frame. This is because the 
fiuid elements can now draw on the potential energy they inherit from starting at the top 
of the bump. 

This suggests that, in the presence of gravity, the results should be substantially different 
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FIG. 4: Same parameters as in Fig. 3, but with initial condition {x^,x'^) = (7r/2,0). For larger 
substrate amplitudes, the system is dominated by 'rimming,' where particles skim a depression 
before moving to the next one, or sometimes undergo long flights. 

if we start elsewhere on the substrate. Figure 4 shows simulations with the same parameters 
as in Fig. 3, but starting at {x^,x'^) — (7r/2,0), some way down the bump. The motion is 
then confined to narrow channels for moderate /q. But for larger substrate amplitudes, the 
system is dominated by 'rimming,' where particles skim a depression before moving to the 
next one, or sometimes undergo long flights. This is a similar situation to basketball (or 
golf), where the ball turns around the hoop (or cup) a while before deciding to go in or 
out. If we take an initial condition at the bottom of the bump, (a;^,a;^) = (tt, 0), then the 
trajectories do not have enough energy to escape the potential well (Fig. 5). 

Finally, in Fig. 6 we illustrate the effect of inclining the substrate at an angle 6 = 7r/8. 
With (f) — —n/2 the trajectories flow 'downhill,' in the negative x'^ direction, modifled by 
the bumps. The system still appears to becomes chaotic for larger /q. Larger bumps 
induce a 'shadow' effect, where they prevent ffuid from ffowing behind them (in particular 
for /o = 0.6). 



V. LYAPUNOV EXPONENTS AND CHAOS 



We now investigate whether the behaviour described in the previous section is chaotic or 
not. The flrst variation of the geodesic Eq. (24) gives 

Sx" + 2Tlp x'^dx^ + d^Tlp x'^x^dx'^ ^ -g F^^^x^ (28) 



10 




FIG. 5: Same parameters as in Fig. 3, but with initial condition = (7r,0). The particles 

begin at the bottom of the potential well, and they do not have the energy to escape. 




FIG. 6: Same parameters as in Fig. 3, but with incline 9 = n/S and 4> = — 7r/2. For larger bumps 
there is a 'shadow' effect, particularly for /o = 0.6. 
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where we have used 9^e'^ — —V^^e^. This equation can be massaged into the geodesic 
deviation equation, [7, 19] 

— + R'p^a x'^x'^Sx^ = (29) 
where D/Dt is defined in (25), and 

R^p-ia — (^i^a(3 ~ ^cX'^p + r^A^a/j " r^;^r^^ (30) 

is the Riemann curvature tensor. For two-dimensional surfaces, the curvature tensor simph- 
fies to 

R^'p-ya — S {S^-y G^a — ^'^ a (31) 

where S = {rt — s'^)/w'^ is the Gaussian curvature and Gjsa is given by Eq. (9). Hence, a 
simphfied form of (28) for surfaces is 

—,Sx'' + 9Ux,x)Sx''-{x,Sx)x'')=0 (32) 

where the inner product is defined by (V, 1^) := GapV°'W^. 

Note that the gravitational term does not enter Eq. (28) directly, though it does 
through (24). For the rest of this discussion we will assume g = 0, since it simplifies 
the discussion considerably. If that is then case, then it is easy to show that 

D D D"^ 

— {x,6x) = {x,—Sx), — (i;,5a;) = (33) 

which means that if we choose the initial Sx" such that {x,Sx) = {x, DSx/ Dt) = 0, 
then {x, Sx) remains zero for all time. With this choice initial condition, the geodesic devi- 
ation equation (32) finally takes the form 

Sx" + 9{x,x)6x'' = 0. (34) 



Now we can ask under what condition the substrate shape will be favourable to chaotic 
geodesies. Since Eq. (34) resembles an oscillator equation, and {x, x) > 0, we see that 
negative Gaussian curvature will favour divergence of trajectories. 

We have not yet solved Eq. (28) for 6x'^, but a comparison of the distance between two 
initially very close trajectories is shown in Fig. 7, for the same parameters as in Fig. 3 
(/o = 0.5). Unsurprisingly, the plot confirms exponential growth, demonstrating at least 
numerically that chaos is indeed present. 



VI. DISCUSSION 

We have shown that the flow of a shallow layer of inviscid, irrotational fluid on a curved 
substrate leads to particle trajectories that follow geodesies in the curved space, subject to 
gravity. We have displayed the range of behaviour that these geodesies can exhibit, from 
regular to chaotic. 

As Fig. 1 shows, the theory is not likely to be valid much beyond the point where char- 
acteristics cross, and viscosity also causes important corrections. Another effect we ignored 
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FIG. 7: The Cartesian distance \6x\ between two trajectories, for the same parameters as in Fig. 3 
with /o = 0.5. The trajectories diverge extremely rapidly, consistent with chaotic behaviour. 

is the possibility that centrifugal forces can cause the fluid to spin out and detach from the 
substrate. [9] Experiments are needed to determine to what extent the chaos observed here 
is reproduced in reality. If chaos is indeed prevalent, then perhaps chaotic advection can be 
exploited in some applications to enhance mixing in shallow layers. 

We made the case in the introduction that chaos in the geodesic equations was a subject 
worthy of study on its own. The emergence of chaotic behaviour as a function of Gaussian 
curvature, as embodied by Eq. (32), should be a rich subject of study, in particular because 
of the simple form this equation takes on a surface. 

We note in closing that a similar study can be made for viscous thin films. [20] However, 
trajectories there are much less prone to chaotic behaviour, because of the diminished role 
of inertia. 
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